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Abstract. Hyperspherical partial wave approach has been applied here in the study 
of double photoionization of the helium atom for equal energy sharing geometry at 20 
eV excess energy. Calculations have been done both in length and velocity gauges and 
are found to agree with each other, with the CCC results and with experiments and 
exhibit some advantages of the corresponding three particle wave function over other 
wave functions in use. 
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1. Introduction 

There has been a significant development in the last one decade in the theoretical study 
of double photoionization (DPI) of the helium atom. This was possible because of 
rapid developments in the experimental side by several groups, extending over several 
countries. Double photoionization of the helium atom is one of the most basic atomic 
processes. Even then, this problem merits further detailed studies as it involves complex 
three body effects and electron correlations that are not yet fully understood. Total 
DPI cross section for helium has been measured by several groups [U El H] and also 
calculated theoretically by many others El HI IE] ■ In general there is good agreement 
between the observed and the calculated cross sections. However, there are also some 
discrepancies. For example experimental results of Samson et al and of Bizau and 
Wuilleumier, although agree within 15 percent of each other, there are systematic 
differences at energies above 1 Rydberg. So also is the case regarding the theories. 
Results of hyperspherical 7£-matrix with semi-classical outgoing waves (H7£M-SOW) of 
Selles et al $ agree favourably with the results of Bizau and Wuilleumier pQ while the 
2SC calculation of Pont and Shakeshaft [3] favours the measurements of Samson et al 
|lj (see ref. As regards the differential cross-sections, understanding of the results 
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has much improved with the availability of results of some elaborate calculations like 
the time-dependent close coupling (TDCC) calculation of Colgan et al [2J, the H.TZM- 
SOW calculation of Malegat et al and Selles et al [HI EI] in addition to the earlier results 
(P33HDH21II3I) and the recent results r|TIl|T51|TB| of CCC theory, the most extensively 
applied theory to the problem to date. For the status of the CCC and some other 
methods the review article by Briggs and Schimdt J7j may be seen. Even then there 
remains the task of improving the theoretical estimates, particularly at low energies and 
for unequal-energy- sharing geometries (for which the results appear comparatively less 
satisfactory) as well as for better understanding of the physics of the DPI problem of 
helium 

In the time-independent frame work the solution of the problem depends basically 
on the accurate computation of the T-matrix element given by 



where ^(ji,^) is the helium ground state wave function, V is the interaction term 
given by 



^i(fi + 7*2) (length form) and Wf if 1^2) is the final channel continuum wave function 
with incoming wave boundary condition for the two outgoing electrons and u>i is the 
incident photon energy. Here e is the photon polarization direction and r\, r*2 are the 
co-ordinates of the two outgoing electrons, the nucleus being at the origin. 

For accurate cross section results one needs accurate wave functions $j and ^ ~\ 
Accurate bound state helium wave functions are easily available. There exists a number 
of such wave functions for the ground state (and low-lying excited states) in analytic 
form of different accuracies, such as simple Hartree Fock type wave function used by 
Maulbetsch and Briggs ^HJ or a Hylleraas type wave function given by Chandrashekhar 
and Herzberg ^H] or by Hart and Herzberg [2U] ■ If needed, one may also readily generate 
an arbitrarily accurate bound state wave function along the line developed by Pekeris 
[2*T] . But for \E'j~' ) there are no such simple accurate wavefunctions in analytical form. 
Most wave functions used in the literature are either incorrect at finite distances or in 
the asymptotic domain, as 3C [22] or 2SC j2H| wave functions. There are now many 
calculations of varied accuracies depending mainly on the use of different final channel 
wave functions. 

Without caring for the explicit form of the wave functions, Huetz and co-workers 
[21] established dependence of cross sections on various angular variables of the outgoing 
electrons and on energy. For the equal energy sharing geometry case it turns out to be 
rather simple in form. On least squares fitting this gives good representation of the 
triple differential cross section(TDCS) results (some authors prefer the name five fold 
differential cross sections (FDCS)which is more appropriate and relevant in view of some 
recent experimental results [25]- However, we will continue to call it TDCS). 



T fi = m >\v\^) 



(i) 



D is the dipole operator and is given D 



V = e-D. 
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For the study of TDCS close to threshold there are the Wannier calculations by 
Faegin 123 ■ These give good representation of the shape of TDCS results at 6 eV 
excess energy but miserably fail at higher energies. 

There are also a number of detailed calculations by Maulbetsch and Briggs fSl EH] 
which used for the final state wave function, the 3C wave function of BBK theory [22] 
and produced moderately accurate cross section results. It is well known that the 3C 
wave function is correct in the asymptotic domain (when all of r 1; r 2 and r 12 are large) 
but it is not accurate enough at finite distances. Similar calculations are reported by 
Pont and Shakeshaft jS] I2S1 EH]- They used screened coulomb (2SC) wave functions 
(for the outgoing electrons) which is supposed to be a better wave function (but not 
asymptotically exact). The results are much better. 

Later Lucey et al jHO] tried various initial state and final state wave functions, 
including the 3C wave function (none accurate enough), tested gauge dependence and 
found much discrepancies in the results. 

In this work, we consider yet another high level computational scheme which is 
capable of yielding reliable cross sections. This is the Hyperspherical partial wave 
approach of one of the authors (Das [21 1321 ESI EH], Das et al (.35.) which is very 
successful in representing the three-particle continuum wave function in the final channel 
very accurately. 

Here we use a 20-term correlated Hylleraas type wave function of the form 

$ i (r 1 ,r 2 ) = Ke-f^ + ^ £ C mi {r x + r 2 ) n \ ri - r 2 p rg (3) 

rii ,n 2 ,n 3 

given by Hart and Herzberg [20J for the helium ground state and use the final channel 
wave function calculated in the hyperspherical partial wave approach as indicated below. 

It may be noted here that our present approach and the H7?.M-SOW approach are 
similar in some respects, although ours is a fully quantal approach whereas the HTZM- 
SOW approach is partly semi classical. From the recent calculations of Selles et al on 
H7£M-SOW approach, it is clear that for unequal energy sharing kinematics one needs 
to consider the asymptotic range parameter (considered below) to have values of 
several thousands a.u. for converged results, whereas for equal energy sharing cases only 
few hundreds suffice. So we restrict our present study to equal energy sharing cases only, 
which is much easier. 



2. Hyperspherical Partial Wave Approach 

In this section we outline the most salient features of this method. For the final 
state ^^f\ which needs more accurate treatment, we use hyperspherical co-coordinates 
R = i/n 2 + r 2 2 , a = atan(r 2 /r 1 ), f\ = r 2 = {0 2 ,(j) 2 ) and u = (a,fi,f 2 ) 

and put P = v / Pi 2 +P2 2 ) "o = atan(p 2 /pi), Pi = (9 Pl ,(p Pl ), p 2 = (O p2 ,(fr P2 ) and 
u = (a ,p\,p 2 ), pi, p 2 being momenta of the two outgoing electrons of energies E 1 
and E 2 and coordinates f[, f* 2 . We expand in hyperspherical harmonics (Das 
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|32| . Lin [HH]) which are functions of the above five angular variables and depend on 
the variables Ei, £ 2 ,n,L,M (collectively called A) which are respectively the angular 
momenta of the two electrons, the order of the Jacobi polynomial and the total angular 
momentum and its projection, in addition to the dependence on S, the total spin. It 
may be noted that L, S, tc (the parity) are conserved here. 

Thus we decompose the symmetrized wave function fy^J as 

= J^E^#("), (4) 
V vr A pi 

following the expansion of the symmetrized plane wave [20] 

[exp(ip~[ ■ f[ + ip 2 - f 2 ) + {-l) s exp(ip 2 ■f[ + ip[- f 2 )]/ (2tt) 3 




'2 



</> s A M. (5) 



Here A = l\ + i 2 + 2 n and p = PR. 

The satisfy an infinite coupled set of equations 



d 2 1 v x {v x + 1 



dp 2 p 2 W v P P 



^ s (p) + E^ k ^(p) = o (6) 



where 



«s* = -mem, (7) 

c = -^-^ + r !— — (8) 

cosa stna \ricosa — r 2 sma\ 

and vx = A + | (note that we use A with two different meanings depending on the 
context). Here the radial waves with L=l, S=0 and tt = — 1 are relevant for the T- 
matrix calculations. 

So here we fix p = (L,S,ir) with L=l, S=0 and n = —1, call it p and consider 
different iV = (£i,£ 2 ,n) and set = f^ in equation (6). Further we omit /x from / 
and write the relevant coupled set of equations (6) as 



N 



2 a s r 



^ + l--^j/i + E^f&-0. (9) 

For our numerical computations we truncate the set to some maximum value N mx of N. 
These N mx equations in N mx variables are needed to be solved from origin to infinity. 
Actually we need construction of N mx independent solutions which vanish at the origin. 
Now for convenience we divide the whole solution domain (0, oo) into three subdomains 
(0, A), (A, Roo) and (-Roo,oo), where A has the value of a few atomic units and Roo 
is a point in the asymptotic domain. Best choices for these may be made by simple 
variations. Results do not depend significantly on these. But for converged results in 
some situations, values of R^, as in H7?.M-SOW calculation are to be thousands of 
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atomic units. Next we proceed for solutions over subdomains. For (Rqo, oo) we have 
simple analytic solutions 



KM = E + d°) 

T P P 

GW = E^ + %^ (id 

£ P P 

where ff® N and are the N-th element of the k-th solution vectors. Obviously 
these give 2N mx independent solution vectors. The coefficients in these expressions 
are determined through recurrence relations (see Das [32] ) in terms of af^ = a^N and 
kfcTV = 0' c i°N = 0' ^kN = a kN, a kN being the N-th element of the k-th eigen vector of 
the charge matrix A = (cynn')- Here we have 9k = p + a^n 2p, a k being the k-th eigen 
value of A. 

Solution over (A, R^) is also very simple. Because of the simple structure of 
equations (9) a Taylors expansion method works nicely. Earlier for the (e, 2e) problem, 
Das also adopted this approach (HHUHl]- But the main difficulty lies in the construction 
of the solution vectors over (0, A). In those calculations on (e, 2e) problem Das used an 
approach as in R-matrix calculations [371 • But very often, this invites pseudo resonance 
type behaviour causing undesirable oscillations in the cross sections. Recently for 
(e, 2e) problem, we applied the finite difference method (a five-point scheme) for 
solutions in the interval (0, A) and thereby get rid of undesirable oscillations. So we 
adopted here the same approach but with a seven-point scheme in place of the five-point 
scheme. 

Thus for the solution in the interval (0, A) we recast equations (9) in terms of R 
instead of p, as 

{d2 +p 2 - — 1) }f N + E f*N> = o. (12) 

N'=l K 



alR 2 R 2 



and solve these equations as a two point boundary value problem by difference equation 
method. At R = 0, the solution vectors are set to zero while at R = A we assign to the k- 
th solution vector the k-th column of the unit matrix. The matrix for the corresponding 
difference equation is a sparse matrix and for its solution special methods are available. 
Here we use biconjugate gradient method [SHI- We find that this method readily works 
and gives converged solutions. 

Now for the difference equations we divide the interval [0, A] into m subintervals 
of length h with mesh points 

= Rq < Ri < R2 < ■ • • < Rk < • • ■ < i? m _i < R m — A 



with R k = Rq + kh and use the following seven-point difference formula: 
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// 1 1 3 

fN( R k) = [ gQ./Ar(-Rfc-3) - 7^16/jv(-Rfc-2) 

3 49 3 3 

+ -^fN(Rk-l) - -^fN(Rk) + ^/aK-Rfc+l) - ^16/Ar(-Rfe+2y 



+ gg./W(-Rfc+3' 



{ 



69 



25200 



(13) 



for A; = 3, 4, • • • , m — 4, m — 3, and for k = 1, 2 and m — 2, m — 1 the formulae 

fN(Ri) = ^[lfN(R ) + Qf N (Ri) -Y h2 fN(Ri) 



and 



h 2 f'^(R 3 ) + 6f N (R 4 + -f N (R 4 ) 



8* 



+ {- 



23 



10080 



h 6 f {viU) (&)}. 



(14) 



h 2 fN(R 4 )+6f N (R 4 ) + -f N (R 5 



8' 



+ {- 

/iv(-Rm-2) = 



23 



10080 
1 r3 



h 2 



fN(Rm-b) + 6/Ar(i? m _ 4 ) — /l f N (R m -4) 



51 



11 



2 

+ {" 



— /l f N (R m -3) + 6/Af(i? m _ 2 ) + ~/iv(-Rm-l 



8' 



23 



10080 



/ l 6/("«)(^ 4 )}. 



(15) 

fN(R m -3) 

(16) 



/jv(-^m-l) 



1 r3 



51 



-jN(Rm-i) + 6/Ar(-R m _ 3 ) — /l f N (R m s) -r/jv(-Rm-2 



11 // 3 

^ /jv(^m-2) + 6/jv(-Rm-l) + giW(-Rm) 



+ {- 



23 



/i 6 / (wi<<) (6)}- 



(17) 



10080 

The quantities on the right hand sides within curly brackets represent the error 
terms. The corresponding difference equations are obtained by substituting these 
expressions the values of second order derivatives from the differential equation (12). 
For continuing these solutions in the domain (A, i?^) we need first order derivatives 
f' N (R) at A. These are computed from the difference formula 



/^J = g^h/iv(i^-4 + 24/ JV (i2 fn _ 2 ) - 128/ JV (i2 m _ 1 ) + 105M/C)] 

+ -=-f N (Ru) + { 



4h A 



(18) 



7 J,,V 105' 

Here too, the quantity within curly brackets represents the error term. The solutions 
thus obtained in (0, A) are then continued over (A, R^) by Taylor's expansion method, 
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as stated earlier, with stabilization after suitable steps j^Hj- The N mx independent 

solution vectors so obtained, are put together to get the solution matrix f . The solution 

(k) 

matrices f sn and f cs are similarly obtained, whose N-kth element are respectively f snN 
and fcsN, given by (10) and (11) respectively. 

Next we introduce the K-matrix through the relation 

fo ■ B = f sn + f cs ■ K (19) 

where B is an unknown constant matrix. The K-matrix is determined from matching 
values and first order derivatives at Roo, where all of fo, f sn and f cs are valid. (It may 
be noted here that there is a slight departure in our definition of K-matrix from the 
usual practice. However, it is symmetric as it should be). 

Finally the physical scattering state with appropriate boundary conditions is taken 

as 

f P h = fo • g (20) 
and also we have 

fph = (fsn + fcs ■ K) ■ C 

= / S n"C + / cs -d (21) 

with 

d = K ■ c. (22) 

Thus the physical state is completely determined once the vector c is determined. Now c 
is determined from the consideration that is asymptotically a (distorted) plane wave 
(representing the two outgoing electrons) plus incoming waves only. So the coefficients 
of the outgoing wave exp(ip) of both ty^J and the symmetrized plane wave (equation 
(5)) must be the same (except for the distorting term exp(iakln 2p)). This requires 

c=[I + iK}- 1 P (23) 

where 

P = -2e*5 X- x $ s *(a;o), (24) 

and X is the matrix comprising of the columns of eigen vectors of the charge matrix A 
and $ s * is given by 

* S *M= ; . (25) 

Finally the DPI triple differential cross section is given by 

d 3 o- = 2ir 2 ap 1 p2 . |2 , , 

after the inclusion of /xo-part of in T/j. 
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3. Results 

In our present calculation we have applied the above hyperspherical partial wave 
approach both in length and velocity gauges. Here we consider equal-energy- sharing 
case, since then the computational problem becomes little simple. For this calculation 
we have chosen A = 5 a.u., Roo = 300 a.u., h = 0.05 a.u. upto A and 0.1 a.u. beyond 
A. We have included 90 coupled channels with n upto 9 and 1%) combinations nearly 
as in ECS calculation [321 for electron - hydrogen ionization collision. We have chosen 
the case of ionization at 20 eV excess energy as it has been widely considered and 
for which there are interesting experimental results jTHj. For the present calculations 
with 90 channels and = 300 a.u., our single differential cross section (SDCS) is 
little above the desired value of about 0.93 Kb/eV at E/2 (E being the excess photon 
energy). So we normalized our TDCS by scaling with a factor 0.8 (which is also the 
factor we use to scale our SDCS to get the desired value of 0.93 Kb/eV at E/2) both in 
the length and in the velocity gauges. The TDCS results thus obtained are presented in 
figure 1. Here we compare our results with the experimental results of Brauning et 
al [TQj and with the theoretical results of the CCC calculation only jTHl Ej, since the 
overall agreement of the CCC results are known to be somewhat better compared to the 
results of other calculations like TDCC, H7?.M-SOW etc. In all the cases the agreement 
between the velocity and length gauge calculations is excellent everywhere. We also did 
calculations in the acceleration gauge and these are indistinguishable from those of the 
velocity gauge. Agreement with the experimental results and with the CCC results are 
also generally good except for some spurious peaks at 9 Pl = 0°. The CCC results appear 
a little better compared to ours. 

For unequal energy sharing our approach also works and we have reasonably good 
results. But for fully converged results, we have to consider larger values for R^. So we 
wish to report such results in future when our study is completed. 

4. Conclusions 

The present calculation, reported here, has only approximately converged. The results 
we have obtained, go to show that the hyperspherical three-particle scattering state 
wave function, used in the present calculation, must be reasonably accurate from small 
distances to the asymptotic region, since the results in all the three gauges are practically 
identical. In contrast, the 3C or other similar wave functions, which are not accurate at 
finite distances, show strong gauge dependence [3U]. We also mention that the present 
calculation is free from any genuine difficulty and does not show any weakness worth 
mentioning. The present approach may easily be applied to DPI with varied types 
of polarization of the incident photons. With judicial choice of the parameter R^ 
and possibly with the availability of better computational facilities, the method may be 
applied from very low energy to high energy cases. At this point if we recall the capability 
of the hyperspherical partial wave approach in representing electron-hydrogen-atom 
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ionization collisions p£H EH at low energies (and also consider situations of very low 
energy cases, with excess energy 1 eV and below for which we had to take about 
4000 - 5000 a.u. and get reliable results for ionization cross sections (121), consider the 
present success, then we may expect the hyperspherical partial wave theory to have a 
very good prospect. 
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Figure Captions 

Figure 1. Triple differential cross sections for photo double ionization of the helium 
atom for equal energy sharing geometry for 20 eV excess energy and for a) 9 Pl = 0°, b) 
9 Pl = 30°, c) 9 Pl = 60°, d) 9 pi = 90°, 9 pi being measured from the photon polarization 
direction. Theory : continuous curve, present calculation in velocity gauge; dashed curve 
present calculation in length gauge; dotted curve, CCC calculation [TU]; Experiment : 
absolute measured values of Brauning et a/[10j. 
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